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Abstract. We investigate proton rapidity correlations for a fireball that fragments due to non-equilibrium 
effects at the phase transition from deconfined to hadronic phase. Such effects include spinodal fragmenta- 
tion in case of first order phase transition at lower collision energies, and cavitation due to sudden rise of the 
bulk viscosity at the crossover probed at RHIC and the LHC. Our study is performed on samples of Monte 
Carlo events. Correlation function in relative rapidity appears to be a sensitive probe of fragmentation. 
We show that resonance decays make the strength of the correlation even stronger. 
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1 Introduction 

Exploration of properties of very hot strongly interacting 
matter is the main goal for which experiments with col- 
liding heavy atomic nuclei at ultrarelativistic energies are 
performed. The suggested phase diagram shows rapid but 
smooth crossover from confined to deconfined matter at 
temperatures around 170 MeV and vanishing baryochem- 
ical potential pj. The smooth crossover becomes a first 
order phase transition at some, presently not precisely 
known, value of the baryochemical potential . To find 
the layout of the phase diagram is one of the main goals 
of present theoretical as well as experimental activity in 
the field. 

The fireball which is created in nuclear collisions ex- 
pands rapidly and so the evolution can proceed out of 
the equilibrium. Particularly, passage of the bulk through 
the phase transition or the crossover can be too fast and 
non-equilibrium scenario could be relevant. In case of the 
first order phase transition this would lead to supercooling 
[3] . Even maximum possible supercooling and subsequent 
spinodal decomposition could happen [Slini[Z]. There, the 
bulk becomes mechanically unstable and fragments into 
droplets of high temperature phase. At RHIC and the 
LHC we likely have the smooth crossover, however. Nev- 
ertheless, in this case calculations show a sharp rise of the 
bulk viscosity [5','5^,'T(11 which might lead to fireball frag- 
mentation [11,12.,13) . as well. Thus the passage through 
the phase transition can in any case lead to fragmentation 
of the fireball into smaller droplets which then evaporate 
hadrons. 

Then, the relevant question is: How does one recognise 
the fragmentation of the bulk? Correlations of hadrons are 
a good candidate for an observable which would be sensi- 
tive to fragmentation. Hadrons originating from the same 



fragment will have correlated velocities. Their velocities 
would be centered around the velocity of the fragment. 
It has been proposed long time ago ITT and reexamined 
more recently [15 that proton correlations are appropri- 
ate observables to look at. Protons are rather abundant 
in nuclear collisions and their velocities are close to the 
velocities of the fragments. This is due to the large proton 
mass which leads to small thermal velocities of protons. 
Pions, in comparison, are much more abundant but also 
very light. Hence, their velocities may differ considerably 
from the velocities of the emitting fragments and their 
correlation function is less suitable to study the effect of 
fragmentation. 



In this paper we study proton correlations in rapidity 
on datasets simulated with a realistic Monte Carlo event 
generator DRAGON [TB]. This is an improvement over 
previous simulations |15j where only a schematic distri- 
bution of droplets has been used. In our model, droplets 
have sizes distributed according to gamma-distribution 
and their positions are random with the only constraint 
that they must not overlap. It is also important that res- 
onance decays are included. We investigate the effect of 
resonances on proton correlations. 



In the following section we define the correlation func- 
tions which have been measured. In Section |3] we give a 
short overview about DRAGON. Results are presented in 
Section m where we also discuss our predictions for FAIR 
and the LHC. Conclusions are given in Section [5] 
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2 The correlation function 

We shall measure the correlation function in kinematic 
variable Ax (e.g. rapidity difference) 



and is normalized to the total number of proton pairs from 
the droplet 
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2, mixed 
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where <p2iAx) is the probability to observe a pair of hadrons 
with the distance Ax. It is averaged over all events of the 
sample. Denominator <?2, mixed (^a:) is the reference distri- 
bution obtained by the mixed events technique. Hence, 
^2, mixed (^a;) is Only determined by the overall distribu- 
tion of hadrons and not affected by any correlations. 
In our study, as Ax we shall use the rapidity difference 



= yi-y2 

and the relative rapidity 
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Note that j/12 is the rapidity that corresponds to relative 
velocity of one particle in the rest frame of the other 
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In the following we re-derive the correlation function 
|15| in Ay for the simple case where protons are emitted 
from droplets distributed according to Gaussian distribu- 
tion. The two proton distribution in rapidity in one event 
is written as a sum of two terms, one corresponding to 
both protons coming from the same droplet, other to each 
proton stemming from a different droplet 



P2(yi,2/2) = X! ^2(2/1 > 2/2) + ^p\{yi)pi{y2) 



(6) 



Here, i,j numerate the droplets and p\ are the single- 
particle distributions for production from one droplet 
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The droplet is thus centred around j/i and its rapidity 
distribution has the width erf, p\{y) is normalized to the 
total number of protons from this droplet Ui 



(8) 



^2(2/11 2/2) dyi dy2 = Vi{vi - 1) ■ 



(10) 



Pl{y)dy = i^i ■ 
The two-proton distribution for one droplet is then 
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The sums in ^ run over all droplets and over all pairs of 
droplets in the event. 

We next introduce two-particle distribution in rapidity 
Terence in one event 



F{Ay) ^ I dyidy2 5{Ay~{yi-y2))p2{yi,y2) 
dyidy2 6{Ay ~ (yi -2/2)) 



This leads to 



XI ^2(2/1, 2/2) + XI ^2(2/1)^2(2/2) 



(11) 



F{Ay) = 



1 1 



exp 



Af_ 

4a2 



X ^i^j exp 



^y - (2/» - Vj) 



(12) 



In order to obtain 'l>2{Ay), this must be now averaged 
over a large sample of events. By doing this we shall also 
sum over many droplets. In the first term on the r.h.s. of 
eq. (fT^ this leads to 



1 
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where the subscript i numbers droplets in one event and 
a stands for different events; {Nd) is the average num- 
ber of droplets per event and {...)m denotes averaging 
over multiplicity distribution for one droplet. Thus (1^(1/ — 
1))m stands for the average number of proton pairs from 
one droplet. Assuming further that droplet centres are 
distributed according to Gaussian distribution with the 
width ^ 
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we can derive that 



^ / ^ N 1 1 
4>2{Ay) = - ^ 
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The denominator of the correlation function is con- 
structed from the mixed pairs sample. Hence, there are 
no pairs where both protons would be emitted from the 
same droplet. The single-particle rapidity distribution is 
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determined by integrating over all possible droplet posi- 
tions 

T{y)= I -=cxp' ^y-^^'^ 
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Note that this distribution is normalized to unity. Thus to 
get ??2.mixod(^y) normalized to average number of pairs 
per event, we calculate 



dyidy2 6{Ay ~ (yi - ya)) t(2/i)t(z/2) 



(17) 



Mixed events are by construction averaged over a large 
number of events and thus we can write directly 



^2,mixcd(^Z/) = -Fraixcd(^y) 

_ l{Nd){iy)Ai{{Nd){:^)M-l) 



exp 
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(18) 

This is almost identical — up to the term involving multi- 
plicity averages — with the second term of the expression 
for <p2{Ay) Uni). In practical evaluation, however, we shall 
normalize the correlation function so that it converges to 
unity for large Ay. This is equivalent to replacing 

{Nd){,y)M {{Nd){iy)M - 1) ^ {NdiNd - l))Ml/ . 
Then we obtain 
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Here we see that t he correlat ion function will have the 
width given by cta/I -|- a'/£,'. Note that ^ is the width 
of the droplet distribution in rapidity, thus the density of 
droplets will be proportional to 



nd 



Based on this, for large {Nd) 
£,{Nd) 
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so the magnitude of the correlation function is inversely 
proportional to the rapidity density of the droplets. This 
is quite natural: In a dense configuration of many small 
droplets — resembling a fog — the system becomes homoge- 
neous and correlation function should be trivial. 

Note that the choice of Gaussian distribution in this 
argumentation is only due to technical simplicity. Quali- 
tatively, the conclusions are more general. 



3 The Monte Carlo event generator 

Monte Carlo events are prepared with DRAGON: Droplet 
and hadron generator for nuclear collisions. Detailed de- 
scription is available in |16| . Here we summarize the main 
features of the model and introduce parameters which are 
varied in subsequent studies. 

Hadrons, including resonances, can be directly pro- 
duced either from the bulk fireball or from the decaying 
droplets. The fraction of those hadrons emitted from drop- 
lets is a parameter that can be set. Resonances decay ac- 
cording to standard two and three-body kinematics with 
constant matrix element. Chain decays are possible. Ex- 
pected dN/dy or the total multiplicity, which is an exper- 
imental observable, is among the input parameters. This 
setting allows for relevant comparisons if the multiplicity 
is kept constant and hadron production is divided between 
bulk and droplet emission in different fractions. 

Directly produced hadrons and the droplets arc formed 
at the usual freeze-out hypersurface of the blast wave 
model. It is characterized by constant longitudinal proper 
time T — \/t' — z' . In this study we only investigate the 
case of central collisions and assume circular transverse 
profile of the fireball with a radius R. The fireball also ex- 
hibits an expansion pattern which is longitudinally boost 
invariant and has tunable gradient in the transverse direc- 
tion. The expansion velocity profile is given as 

u^{x) — u^{r\, r, (j)) 

— (cosh 77 cosh?7f(r), cos(j) sinhr7t(r), 

sine/) sinh77f(r), sinhr/ cosh 774 (r)) (22) 

(23) 



?7t(r) = \/2po-^ ■ 



where we use the polar coordinates r and (j) in the trans- 
verse plane and the space-time rapidity 

1 , t + z 

•n — — in . 

' 2 t- z 

These coordinates are taken in the centre of mass system 
of the collision. 

Droplets are randomly placed within the fireball. Their 
volumes follow the Gamma distribution 



V2{V) 



V 
6^' 



-V/b 



(24) 



where the volume parameter b determines the average vol- 
ume 

The masses of the droplets are proportional to their vol- 
umes while the energy density is taken with the value 
0.7 GeV.fm"'^. Their velocities are given by u^^{x) at the 
position where they are produced. Droplets decay expo- 
nentially in time with the mean lifetime being equal to the 
radius of the droplet. 

Hadrons, regardless if they come from a droplet or from 
bulk, are emitted from a thermal source with kinetic tem- 
perature Tkin- In case of production from bulk, the thermal 
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distribution is boosted with the local expansion velocity at 
the point of hadron emission. In case of production from 
a droplet the thermal distribution is assumed in the rest 
frame of the droplet. 

Chemical composition of the produced hadrons is de- 
termined by the assumption of chemical equilibrium with 
temperature Tch- Chemical potentials correspond to con- 
served quantum numbers with baryonic chemical potential 
Hb being a free parameter and strangeness chemical po- 
tential being set by the requirement of total strangeness 
neutrality. 

Resonances with masses up to 1.5 GeV/c^ (mesons) 
and 2 GeV/c^ (baryons) according to the particle data 
book are included. 
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4 Results 

Correlation functions are evaluated from samples of artifi- 
cial events. Numerators are built as histograms including 
pairs of protons from the same event according to their 
relative rapidity. We investigate rapidity difference in 1 
dimension, eq.([21), and the relative rapidity in three di- 
mensions, eq.(j31). Denominator of the correlation function 
is constructed from pairs of randomly chosen protons from 
different events. The constructed correlation functions are 
normalized so that they approach unity for large rapidity 
differences. The obtained data are fitted by Gaussian form 
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where x is replaced by either yi2 or Ay and A, a are fit 
parameters. 

Most detailed studies are performed on samples of events 
generated with RHIC parameter settings. Chemical com- 
position is characterized by Tch = 170 MeV and /is = 
46 MeV. Total hadron rapidity density is dN/dy = 1000. 
(Note that this includes also neutral hadrons.) Kinetic 
temperature, if not stated otherwise, is put to 150 MeV. 
There is transverse fiow with po — 0.7. Rapidity distribu- 
tion is uniform, but only hadrons from the interval [—1, 1] 
are taken into the analysis. For each parameter settings a 
sample of 10,000 events is generated. 

We first consider the case of only partial production 
of hadrons (protons) from the droplets. In Fig. [T] we keep 
the average size of droplets constant by fixing 6 = 50 fm 
and change the fraction of hadrons being emitted by the 
droplets. The other part comes from the bulk fireball while 
the multiplicity is kept constant. Clearly, at vanishing d 
the correlation function is trivial. In both cases, for cor- 
relation functions in ID variable Ay and 3D variable j/12, 
the width is independent of the fraction d while the cor- 
relation strength A increases as more hadrons come from 
the droplets. We also see that the correlation is about a 
factor of 2 stronger if plotted as a function of yi2 (3D cor- 
relation). Therefore, in the following we shall only study 
correlation function in j/12. 

The effect of varying the size of droplets is investigated 
in Fig. [31 Here, one half of all hadrons comes from droplets. 



Fig. 1. Correlation functions in yi2 (upper panel) and Ay 
(lower panel) for RHIC parameter settings with b = 50 fm^ 
and varying fraction of droplet-emitted hadrons: d = (black 
dash-dot-dotted curves), d — 0.25 (green solid), d = 0.5 (blue 
dash-dotted) , and d = 1 (red dotted) . 



Since large droplets correlate more protons, the strength 
of the correlation function is increased as b grows. The 
correlation goes away when the droplets are so small that 
not enough energy is available for a proton pair creation. 
At 6 = 10 fm^ the average mass of a droplet is 3.5 GeV 
which may be enough for two protons to be created, thus 
we still see correlation at the level of 0.01. At b = 5fm'^ 
this decreases to 1.75 GeV which is insufficient for a proton 
pair and the correlation function would be trivial. 

We observe that the correlation strength grows both 
with bigger droplets and/or with larger droplet fraction. 
The dependence of A on 6 and d is studied in Fig. [31 Here, 
we summarize the values which were obtained from fits to 
correlation functions. We conclude that it is not possible to 
determine the droplet size or the droplet fraction just from 
measuring the correlation strength A. The value of A can 
at most help to establish a range of possible combinations 
of b and d. If, however, one of them can be determined by 
some other means, then rapidity correlations can be used 
for the measurement of the other. 
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Fig. 2. Correlation functions in yi2 for RHIC parameter 
settings with one half of hadrons coming from droplets, i.e. 
d = 0.5. Different curves correspond to different sizes of the 
droplets: b = 10 fm^ (black dash-dot-dotted curve), b = 20 fm^ 
(green solid), b — 50 fm^ (blue dash-dotted), and b — 100 fm^ 
(red dotted). 



Fig. 4. Correlation functions for RHIC parameters with all 
hadrons emitted from droplets and & = 50 fm"^ . Different curves 
correspond to different kinetic freeze-out temperatures: T = 
170 MeV (black dash-dot-dotted curve), T = 150 MeV (green 
solid), T = 130 MeV (blue dash-dotted), T = 110 MeV (red 
dotted). 
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Fig. 3. Contour plot shows the dependence of the amplitude 
A of the correlation function calculated for RHIC as a func- 
tion of droplet size parameter b and fraction of particles from 
droplets d. The step between the contours is 0.05, labels show 
the corresponding value of the contours. 



So far we have only investigated the variation of the 
correlation function with changing properties of the drop- 
lets. In Fig. m we study the effect of temperature variation. 
Average droplet volume is given by 6 = 50 fm'^ and all 
hadrons are emitted from droplets. We observe weakening 
of the correlation as the temperature increases. The inter- 
pretation is at hand: the effect of the temperature is the 
smearing of hadron velocity so that it does not exactly 
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Fig. 5. Correlation functions in yi2 for standard RHIC pa- 
rameter settings with all particles emitted from droplets and 
b = 50 fm^ Different curves show: full correlation function with 
all resonance decays included (red dotted curve), result from 
a simulation with no resonance production included (green 
solid), result from a simulation with resonances included but 
protons coming from decays of Delta resonances were not taken 
into analysis (blue dash-dotted). 



follow that of the droplet. Therefore, lower temperature 
leads to increased correlations. Formally, this effect can 
be also seen in eg. dT^ (although that equation was de- 
rived for Ay). Higher temperature translates into larger 
(T. Through the pre-factor 1/a (at large this makes the 
correlation weaker. 

Our model also allows for dedicated investigation of 
the effect of resonance decays on the observed correla- 
tions. This is shown in Fig. [S] We see that the resonances 
increase the intercept of the correlation function, i.e., they 
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strengthen the correlation. This is a result of two effects: in 
resonance decay, daughter particles obtain some momen- 
tum due to the mass difference between the resonance and 
its decay products. This would lead to more smearing of 
the final hadron momenta and lower correlation. On the 
other hand, resonances are heavier and therefore closer to 
the velocity of the emitting droplet than lighter hadrons. 
This increases the correlation. Eventually, the second ef- 
fect wins and the correlation function is increased. Note 
that the smearing of resonance momentum in production 
from the droplet gets weaker if the temperature is low- 
ered. Since in our simulation we calculate with rather high 
kinetic freeze-out temperature of 150 MeV, it can be ex- 
pected that this effect will always dominate the influence 
on the correlation function. 

In Fig. [S] we also see that the dominant contribution 
within the resonance decays that produce protons comes 
from the deltas. 

We next turn to predictions of the proton rapidity cor- 
relations at other collision energies. First is that of few 
GeV per nucleon. This is the domain that should be stud- 
ied with the help of the planned programs at low energy 
RHIC run, SPS energy scan, FAIR in Darmstadt, and 
NICA at the JINR Dubna. This energy domain seems 
important since the critical point of the QCD phase di- 
agram is suspected to be here. If so, these accelerators 
might bring us into the region of the phase diagram with 
first order phase transition where spinodal decomposition 
[SIE] and subsequent emission of hadrons from droplets 
can happen. 

For FAIR energy we choose Gaussian profile of the 
fireball in the space-time rapidity 77, with the width Arj = 
0.7. Chemical composition is taken from the analysis of 
Pb-I-Pb collisions at CERN SPS at projectile energy 40 A 
GeV [T7] : chemical freeze-out temperature is 140 MeV and 
the baryochemical potential is 428 MeV. Kinetic freeze- 
out temperature agrees with that of chemical freeze-out. 
Again, there is transverse flow with po — 0.4. The total 
multiplicity is 1250. As a result, on average there are 135 
protons per event. Results on the correlation functions are 
very similar to those obtained for RHIC energy and are 
therefore not shown separately. 

We also make simulations for a flreball as expected in 
Pb+Pb collisions at the LHC. Due to near nuclear trans- 
parency the baryochemical potential is taken at 1 MeV 
while the chemical freeze-out temperature is 170 MeV. Ki- 
netic freeze out temperature is set to 150 MeV and there 
is slightly higher transverse expansion gradient po = 0.8. 
Multiphcity density dN/dy = 2000. 

Since results at all investigated collision energies are 
qualitatively very similar, we only show in Fig.[S]a compar- 
ison of the correlation functions at different collision ener- 
gies. The amplitude of the correlation function is highest 
for FAIR. This is due to the regime with lowest average 
dN/dy where clusterization in momentum space is best 
visible. Somewhat surprising, the FAIR correlation func- 
tion is also the widest. We expected just the opposite since 
the kinetic freeze-out energy at FAIR was lower than in 
the other two cases. The larger width of the correlation 
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Fig. 6. Comparison of correlation functions in yi2 for different 
collision energies. For parameter sets see text. All hadrons are 
emitted from droplets with b = 50 fm'^. Different curves show: 
FAIR (red dotted curve), RHIC (blue dash-dotted), RHIC 
without resonances in the simulation (green solid) , LHC (black 
dash-dot-dotted) . 

function is caused by the Gaussian profile of the fireball 
in space time rapidity in comparison with uniform profiles 
at higher energies. Again, we can understand the effect to 
some extent from eg. (1191) : a source with narrow rapidity 
distribution will have smaller ^. This broadens the peak of 
the Gaussian correlation function. Paradoxically, fireball 
that is very narrow in space-time rapidity will have wide 
correlation function. 



5 Conclusions 

Correlation functions in relative rapidity appear to be a 
powerful tool for the study of clustering of hadron pro- 
duction due to fragmentation. Note that in calculating 
the correlation functions we did not include the effect of 
final state interaction (strong and Coulomb) between the 
protons |14) , which governs the correlation function in ab- 
sence of other effects. This method was also followed in 
|15) . Hence, our results cannot be directly compared to 
experimental data and should be regarded as an estimate 
of the size of the effect. It is large indeed, if the droplets 
have average volume of few tens of fm'^. With the average 
size reaching down to the level of a few fm'^ the infiuence 
on the correlation function is on the per cent level. This is 
the size of the droplets which may be expected at RHIC, 
as indicated by data from PHOBOS. At the LHC, the 
droplets, if any, are likely to be smaller [XT'. 

This study assumed that after emission from the droplets, 
hadrons do not rescatter. Rescattering would weaken the 
correlation and might even wash it away. 

Also, note that we did not study correlation of three 
and more protons 115] . Such correlations demand much 
better statistics and are small if the average volume of 
droplets drops considerably below 10 fm^. For such small 
droplets it becomes unlikely that more than two protons 
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will be emitted from one droplet. Note that we do not 
expect at RHIC droplets much bigger than this (if any). 
Thus two-particle correlation appears to be the most suit- 
able tool at these collision energies. 
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